home *** CD-ROM | disk | FTP | other *** search
/ Ham Radio 2000 #1 / Ham Radio 2000.iso / ham2000 / misc / radiosky / azel.bas next >
Encoding:
BASIC Source File  |  1996-05-15  |  5.1 KB  |  126 lines
  1. 1 'The following simple program, written in IBM BASIC, will perform the
  2. 2 'calculations to convert HA, decl to and from az-el, including the
  3. 3 'computation of LST. Although this program prints out the results
  4. 4 'to limited precision, the calculations are precise to a small
  5. 5 'fraction of a degree, and are valid at least throughout this and
  6. 6 'the next century. SKYCOORD.BAS is a complete working program, but
  7. 7 'the reader may wish to use it as a template for his or her own program
  8. 8 'in any favorite programming language, perhaps adding a more sophisticated
  9. 9 'user interface.
  10.  
  11.  
  12.  
  13. 10 REM Convert RA, Dec to Az, El and vice versa.
  14. 20 REM    Put together by Darrel Emerson, aa7fv, July 1995.
  15. 25 REM    From various sources, especially "Astronomy with your PC"
  16. 26 REM    by Peter Duffet-Smith, Cambridge University Press. 
  17. 30 DEFDBL A-Z
  18. 40 pi = 4# * ATN(1#)
  19. 50 INPUT "Day, Month, Year (e.g. 17,6,1995)"; DY, MN, yr
  20. 70 INPUT "Observer's longitude (decimal degrees, W is negative)"; lo
  21. 80 gl = lo * pi / 180#
  22. 90 INPUT "Observer's latitude (decimal degrees, S is negative)"; la
  23. 100 tl = la * pi / 180#
  24. 120 PRINT
  25. 130 PRINT "RA-Dec to Az-El (1)": PRINT "Az-El to RA-Dec (2)"
  26. 140 INPUT " or finish.     (3)   "; ir
  27. 150 IF ir < 1 OR ir > 2 THEN PRINT "End conversion program": STOP
  28. 160 PRINT : INPUT "UT hours (decimal hours: e.g. 23:45 = 23.75)"; ut
  29. 170 GOSUB 2100: REM calculate LST hours
  30. 180 PRINT "LST"; USING "####.####"; sl; : PRINT " hours"
  31. 190 IF ir = 1 THEN GOSUB 300
  32. 200 IF ir = 2 THEN GOSUB 500
  33. 210 GOTO 120
  34.  
  35. 300 REM Convert RA-Dec to Az-El
  36. 310 INPUT "Right Ascension, decimal hours (e.g. 20h54m=20.9)"; rh
  37. 320 INPUT "Declination, decimal degrees"; dd
  38. 330 de = dd * pi / 180#
  39. 340 ha = (sl - rh) * 15# * pi / 180#
  40. 350 GOSUB 3000
  41. 360 PRINT "Azimuth:  "; USING "####.#"; az * 180# / pi; : PRINT " degrees"
  42. 370 PRINT "Elevation:"; USING "####.#"; el * 180# / pi; : PRINT " degrees"
  43. 380 RETURN
  44.  
  45. 500 REM Convert Az-El to RA-Dec
  46. 510 INPUT "Azimuth, decimal degrees"; ad
  47. 520 INPUT "Elevation, decimal degrees"; ed
  48. 530 az = ad * pi / 180#
  49. 540 el = ed * pi / 180#
  50. 550 GOSUB 4000
  51. 560 hh = ha * 180# / pi / 15#: REM convert hour angle into hours
  52. 570 REM print"Hour Angle:";using "###.####";hh;:print " hours"
  53. 580 ra = sl - hh
  54. 590 IF ra < 0# THEN ra = ra + 24#
  55. 600 IF ra > 24# THEN ra = ra - 24#
  56. 610 PRINT "Right Ascension: "; USING "###.##"; ra; : PRINT " hours"
  57. 620 PRINT "Declination:     "; USING "###.#"; de * 180 / pi; : PRINT " degrees"
  58. 630 RETURN
  59.  
  60. 1000 REM subroutine atan2, using sine s and cosine c.
  61. 1005 REM All angles radians.
  62. 1010 IF c = 0 AND s = 0 THEN atan2 = 0: RETURN
  63. 1020 IF c = 0 AND s > 0 THEN atan2 = pi / 2#: RETURN
  64. 1030 IF c = 0 AND s < 0 THEN atan2 = -pi / 2#: RETURN
  65. 1040 atan2 = ATN(s / c): IF c < 0 THEN atan2 = atan2 + pi
  66. 1050 IF atan2 > pi THEN atan2 = atan2 - 2# * pi
  67. 1060 REM atan2 is returned between +/-pi
  68. 1070 RETURN
  69.  
  70. 1100 REM subroutine to calculate Julian days since the start of the
  71. 1105 REM century. (i.e. since 1900 Jan 0.5)
  72. 1110 REM TAKES YR (E.G. 1995), MN (E.G. 6), DY (E.G. 17) & gives DJ
  73. 1120 M1 = MN: Y1 = yr: B0 = 0
  74. 1130 IF MN < 3 THEN M1 = MN + 12: Y1 = Y1 - 1
  75. 1140 A0 = INT(Y1 / 100): B0 = 2 - A0 + INT(A0 / 4)
  76. 1150 C0 = INT(365.25# * Y1) - 694025!
  77. 1160 D0 = INT(30.6001# * (M1 + 1)): dj = B0 + C0 + D0 + DY - .5#
  78. 1170 RETURN
  79.  
  80. 2100 REM Find GST (=sg) and LST (=sl) from UT,
  81. 2105 REM using UT, YR, MN, DY, GL (=longitude, W neg)
  82. 2110 GOSUB 1100: REM get DJ, days since 1900
  83. 2120 d = INT(dj - .5#) + .5#: t = (d / 36525#) - 1
  84. 2130 r1 = 6.697374558# + (2400# * (t - ((yr - 2000#) / 100#)))
  85. 2140 r0 = t * (.0513366# + t * (.00002586222# - t * .000000001722#))
  86. 2150 t0 = (r0 + r1) - 24# * INT((r0 + r1) / 24#)
  87. 2160 sg = (ut * 1.002737908#) + t0
  88. 2170 sg = sg - 24# * INT(sg / 24#)
  89. 2180 sl = sg + (gl * 180# / pi / 15#):
  90. 2190 sl = sl - 24# * INT(sl / 24#)
  91. 2200 REM sg is Grenwich Sidereal Time, sl is Local Sidereal Time
  92. 2210 RETURN
  93.  
  94. 3000 REM Find Azimuth AZ and Elevation EL.
  95. 3010 REM from Hour Angle HA, Declination DE and Terrestrial Latitude TL.
  96. 3020 REM All angles radians. 
  97. 3030 c1 = -COS(de) * COS(ha) * SIN(tl) + SIN(de) * COS(tl)
  98. 3040 c2 = -COS(de) * SIN(ha)
  99. 3050 c3 = COS(de) * COS(ha) * COS(tl) + SIN(de) * SIN(tl)
  100. 3060 s = c2: c = c1: GOSUB 1000: az = atan2
  101. 3070 s = c3
  102. 3080 IF ABS(SIN(az)) > .7 THEN c = c2 / SIN(az): GOTO 3110
  103. 3090   c = c1 / COS(az)
  104. 3100   REM The "0.7" is a criterion for minimizing numerical rounding error
  105. 3110 GOSUB 1000: el = atan2
  106. 3120 IF az < 0 THEN az = az + 2# * pi
  107. 3130 REM "az" now between 0 and 2*pi, with Az=0 north
  108. 3140 RETURN
  109.  
  110. 4000 REM Find Hour Angle HA, and Declination DE,
  111. 4010 REM from Azimuth AZ and Elevation EL, with Latitude TL.
  112. 4020 REM All angles radians.
  113. 4030 c1 = -COS(el) * COS(az) * SIN(tl) + SIN(el) * COS(tl)
  114. 4040 c2 = -COS(el) * SIN(az)
  115. 4050 c3 = COS(el) * COS(az) * COS(tl) + SIN(el) * SIN(tl)
  116. 4060 s = c2: c = c1: GOSUB 1000: ha = atan2
  117. 4070 s = c3
  118. 4080 IF ABS(SIN(ha)) > .7 THEN c = c2 / SIN(ha): GOTO 4110
  119. 4090   c = c1 / COS(ha)
  120. 4100   REM  The "0.7" is a criterion for minimzing numerical rounding error
  121. 4110 GOSUB 1000: de = atan2
  122. 4120 IF ha < 0 THEN ha = ha + 2# * pi
  123. 4130 REM ha now between 0 and 2*pi
  124. 4140 RETURN
  125.  
  126.